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AN ENERGY DECAYING SCHEME FOR NONLINEAR DYNAMICS OF SHELLS 

CARLO 1.. BOTTASSO\ OLTVTER A. BAUCH AU t 1 AND JOU- YOUNG CHOI* 


Abstract. A novel integration scheme for nonlinear dynamics of geometrically exact shells is developed 
based on the inextcnsible director assumption. The new algorithm is designed so as to imply the strict decay 
of the system total mechanical energy at each time step, and consequently unconditional stability is achieved 
in the nonlinear regime. Furthermore, the scheme features tunable high frequency numerical damping and 
it is therefore stiffly accurate. The method is tested for a finite element spatial formulation of shells based 
on mixed interpolations of strain tensorial components and on a two-parameter representation of director 
rotations. The robustness of the scheme is illustrated with the help of numerical examples. 

Key words, geometrically exact shell, time integration, energy decaying scheme 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. The formulation of integration algorithms for nonlinear dynamics of geometrically 
exact shells is the focus of this work. The partial differential equations governing this class of problems are 
known to present a rich mathematical structure. In particular, the resulting models are Hamiltonian systems 
characterized by a symplectic nature and associated with conservation laws that stem from symmetries of 
the Hamiltonian. The linear and angular momentum as well as the total mechanical energy are conserved 
for free motions of such systems. 

The understanding of the geometric characteristics of the governing equations has been historically 
confined to the fields of analytical mechanics and pure mathematics. Surprisingly, this knowledge has been 
seldom used for the development of numerical methods. Indeed, the study of new integration algorithms has 
been traditionally preoccupied with the development of methods applicable to vast classes of problems, for 
example the class of differential/algebraic equations, or hyperbolic conservation laws. Consequently, classical 
methods rarely preserve the underlying structure of the problem being solved, and hence, such structure is 
lost in the numerical solution. 

This approach also limits the possible theoretical analyses of the schemes, which are, more often than 
not, confined to linear or model cases. For instance, it is customary to characterize integration schemes 
for structural dynamics by studying their behavior when applied to a linear oscillator. This approach is 
clearly not adequate when dealing with highly nonlinear problems such as the dynamics of geometrically 
exact shells. 

A new approach to the design of integration algorithms attempts to bridge the divide between theoretical 
and numerical mechanics. Under this new paradigm, numerical schemes are “backward-engineered” to 
preserve some important qualitative features of the governing equations. Fittingly, this approach is now 
called geometric integration in the mathematical community [18]. Attempts at designing “geometry- aware” 
algorithms for structural dynamics problems can be traced back to the work of Sinio and co-workers who 
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analyzed the problems of rigid body dynamics [26], nonlinear elasto-dynamics [23]. geometrically exact 
shells [24], and geometrically exact beams [25]. In all cases, the idea was to design algorithms that ensure 
the discrete preservation of the total mechanical energy and of the linear and angular momenta of the system. 

When integrating linear and nonlinear finite element models, the implications of the discrete equations 
stiffness must be carefully considered. Indeed, high frequencies are an artifact of the spatial discretization 
process and do not reflect the high frequencies of the original infinite dimensional problem. The need for 
high frequency numerical dissipation has been recognized in the past for linear problems [20]. When dealing 
with complex nonlinear systems, numerical dissipation becomes indispensable. Indeed, nonlinearities provide 
a mechanism for transferring energy from low to high frequency modes. Consequently, numerical solutions 
feature violent oscillations of a purely numerical origin that will eventually play havoc with the convergence 
characteristics of the nonlinear equation solver. 

Among the various geometric characteristics of shell equations, energy preservation appears to be the 
most important for the development of robust time integration schemes. In fact, strict energy preservation 
at the discrete level leads to unconditional stability in the nonlinear regime , whereas the classical approach 
based on the analysis of the spectral radius leads to unconditional stability in the linear regime only. An 
energy preserving (EP) scheme for geometrically exact shell is developed in this paper. In addition, the 
scheme also preserves both linear and angular momenta of the system at the discrete level. Unfortunately, 
preservation of energy and high frequency dissipation cannot coexist, unless energy is transferred from high 
to low frequency modes, a transfer that has no physical basis. To solve this problem, a family of energy 
decaying (ED) schemes that imply a controllable energy decay within each time step is proposed in this work. 
In geometric terms, this means that the evolution of the system is not confined to the level set of constant 
energy, but is allowed to drift away from it in a monotonic and controllable manner. Since the energy remains 
bounded at all times, the scheme is unconditionally stable for nonlinear systems. Furthermore, it can be 
shown that the energy dissipation mechanism of the algorithm is the result of the removal of the higher 
frequencies from the computed response. 

In related papers, various energy preserving and decaying geometric integrators were developed for rigid 
bodies and geometrically exact beams [10, 11, 14, 6, 7], and nonlinear elastodynamics [9]. The concepts 
were extended to multibody systems featuring nonlinear holonomic constraints [3, 15, 12, 16, 13]; non- 
holonomic and unilateral constraints were treated in [5, 4]. The integration of the present shell model in a 
general finite element based multibody framework is discussed in [8]. The proposed scheme is independent 
of the choice of spatial discretization applied to the governing partial differential equations. In the present 
implementation, the finite element method is used, and the mixed interpolation of tensorial components [1, 
2, 17] is implemented to avoid the shear locking problem. The orientation of unit shell directors is described 
by a special family of two-parameter rotations. 

The paper is laid out as follows. The classical equations of motion for geometrically exact shells based 
on inextensible unit directors are presented in section 2. Next, an EP scheme is developed in section 3. 
Section 4 then presents an ED algorithm with tunable high frequency dissipation that is constructed from 
the EP scheme. Finally, numerical examples are presented in section 5 to demonstrate the efficiency and 
robustness of the proposed scheme. A discussion section concludes the paper. 

2. Formulation of the Equations of Motion. 

2.1. Kinematics of the Shell Problem. Consider a shell of thickness h and reference surface area fi, 
as depicted in fig. 2.1. An inertial frame of reference S consisting of three mutually orthogonal unit vectors 
ii, i 2 , i 3 is used. Let r 0 be the position vector of an arbitrary point on the reference surface of the shell, and 
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let C be the material coordinate along n, the normal to the reference surface. The position vector r of an 
arbitrary point on the shell in its reference configuration is then 

r(£\£ 2 , o - EoK'^ 3 ) + C *(C,£ 2 ), (2.1) 


where £ ] a.nd £ 2 are the material coordinates used to represent the shell reference surface. The coordinates 
, £ 2 and ( form a set of curvilinear coordinates that are a natural choice to represent the shell geometry. 
The coordinates £ ] and £ 2 are assumed to be lines of curvatures of the shell reference surface. The base 
vectors are then 


9 = 



9 r 



dr 


r r 

-* 1 ’ -> 2 ’ dc. 

— 

(1 - flj, (1 - -£-) n 2 . n 

n 2 


( 2 . 2 ) 


where 7?i and R 2 are the principal radii of curvature, a n = r# , and the notation (-) |Q is used to denote a 
derivative with respect to It is convenient to introduce a set of three mutually orthogonal unit vectors 
at the shell reference surface (i.e. at £ = 0) 







= ZL, 


(2.3) 


where n aa = a a ■ a a . 

Two fundamental assumptions will be made concerning the deformation of the shell, i.e. the material 
line initially normal to the reference surface of the shell remains a straight line and suffers no extension. 
This is the classical inextensible director model. With these assumptions, the position vector of a material 
point of the shell writes 


me , <e 2 , o - zr>(£ 1 1 , ^ 2 ) + Hie ,e)+cE,(t\ eh 


(2.4) 
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whore «(£' ,£ 2 ) is the reference surface displacement vector. In the deformed configuration, the base vectors 
at. the shell reference surface are 


G = [G„ Ga, &] = 



R_ 2 i 


dR' 

3C. ' 


(2.5) 


Introducing the position vector, eq. (2.4), then yields 

d= [.a. 3 l _ 

V5TT’ y/a£ 


. Ga 


— E + ( //, 


( 2 . 6 ) 


where 


b = i 2 , &] = 


£, + 


v^TT 


£9 + 


u 2 


H = 


E, 


E 


.3,2 


LV^TT’ \/ a 22 


,0 


(2.7) 


Note that E^f 1 ,*; 2 ) is a unit vector, whereas E^ and E 2 are not unit vectors, nor are they orthogonal to 
E 3 , as axial and transverse shearing strains develop din ing deformation. 


2.2. Equations of Motion. The Green- La grange strain tensor e is defined as 

e = \ (G t G - g T g). (2.8) 

The strain tensor e is defined in the curvilinear coordinate system defined by coordinates £ l ,£ 2 and 
However, it is more convenient to work with the strain tensor e defined in the locally rectangular system 
defined by triad e } , e 2 , see e QS. (2.3). For shallow shells ( i.e . £/i?i <§C 1 and C /^2 ^ 1) undergoing large 
displacements and rotations but small strains (all strain components are assumed to be small compared to 
unity), the strain -displacement relationships can be written as 

e = l[E T E- ■I + C(E t H+H t E + k)], (2.9) 

where 


1/E, 0 0 

0 1 /R 2 0 

0 0 0 


(2.10) 


It is clear that the strains can be expressed in terms of five parameters: the three components of the 
displacement field u (through E, and E 2 ) and the two parameters defining the orientation of the unit 
director E 3 . Virtual changes in the strain energy of the structure are given by 

SV = [ f SV dCdfl = [ f fie * t dCdH, (2.11) 

Jo Jh Jo J h 

where SV is the virtual strain energy density, and r the second Piola-Kirchhoff stress tensor. Introducing 
the strains, eq. (2.9), and taking into account the symmetry of the stress tensor then yields 


SV = 6E-{E + C H)t + SH - (Er. 


( 2 . 12 ) 


The existence of a strain energy density function V is postulated here, hence the constitutive laws are of the 
form r — dV/de. 

The velocity vector of material point P of the shell is obtained by differentiating the position vector, 
eq. (2.4), with respect to time, to find v — u + (E 3 . The kinetic energy of the system is now 

K= f f EdCcin = l [ [ pv-v dCdn, (2.13) 

Jo Jh 2 Jq J h 


•1 


(2.14) 


where K is the kinetic energy density. Introducing the velocity vector then yields 

K ~ 2 P ^ + ^ 3 ) ‘ (ll + (£ 3 )- 

Hamilton T s Principle can now be expressed as 


f tf [ [ (sit - w 

Jti Jn Jh 


) d(dQd t 


r f r 

Jti Jq Jh L 


p(Su + C SE 3 ) ■ (fi + C£ 3 ) + 8E • (E + (H)t + 8H • C^r] d(dQd* 


0 . (2.15) 


Integrating through the thickness of the shell, we get 


[ ' J {% [h-(Ki, i + K- 2 . 2 )] + SE ■ [g f - (A£, t , + M 2 , 2 ) + 2 ^ 3 ] } dndf 


= 0. 


(2.16) 


In this expression, h = mw + s*E$i and g — s*u + I*E$ are the linear and angular momentum vectors of 
the shell, respectively; the mass coefficients are defined as m = J h p d£, s* = f h pQ d(, /* = f h pC d(\ 
The in-plane forces are N_ a — ( EAT * + HM^)/ x /a nn , the out-of-plane forces N$ = E N 3 , and the bending 
moments M n = The convected forces are A r * = [JV* , 7V.^ , 7V£] = J h r d£, and the converted 

bending moments M* = [M * , M * , M 3 ] = f h r( d£. 

The equations of motion of shells could be derived from this principle by expressing the variations 8E 3 
in terms of two components of virtual rotation. 

3. Energy Preserving Scheme. Discrete equations of motion that imply discrete conservation laws 
for the total mechanical energy, linear momentum and angular momentum of the system will now be de- 
veloped. Times £* and t j denote the initial and final times for a time step, respectively, and the subscripts 
(-)i and (•)/ indicate quantities at ti and tj, respectively. Furthermore, the subscript (*) m is used to denote 
mid-point average quantities defined as 

(Om = \ [(■)/ + (•)*] ■ (3-1) 

The following matrix identity will be used extensively 


AjB f -AjBi = {A s - Ai) r B m + A T m {B } - B t ). 


(3.2) 


Hamilton’s Principle, eq. (2.15), is now approximated in time in the the following manner 



fn Jh 


P [(«/ - U t ) + C(i? 3 / - # 3 i)] 


Lf ~ Mi Jk.y ~ &3i 


At At 

+{Ej - Ei) ■ (E m + < H m )r n + ( Hf - Hi) - (E m T n i dC<ifi - 0. (3.3) 


The change in strain components from ti to tj is evaluated with the help of identity (3.2) to find 

ej - Pi = \ [{E f - Ei) T {E m + C H m ) + [E m + C H m ) T (E f - E.) 

+ (H f - Hi) T CE m + C Ej n {Hj - Hi)] . 


(3.4) 
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Over one time step, the strain components can be approximated as e(r}) — e m -f r}(ej - e,)/ 2 , where 77 = 
2(f, - t TT ,)/At is the non-dimensional time. If the strain energy density function V is viewed as a function of 
the scalar variable 77 , the mean value theorem then implies the existence of a fj G [ — 1 , 1 ] such that 


l/ = Vi + 


at; - , . 

—2 = Vj + r n - (e.f - c,). 
n dr > 


(3.5) 


This relationship defines the average second Piola-Kirchhoff stress tensor, r a — dV/de \ Combining this 
result with eq. (3.4) then leads to 


(E f - Ei) • (E m + CH m )r a + (F/ - F,) ■ <F m r n - (e f - e,) ■ r rt = V f - V u (3-6) 

where the symmetry of the stress tensor was taken into account. For linear constitutive laws of the form 
t = C* e, where C* is the stiffness matrix, the average stress tensor simply becomes r a = C* e m . 


The following configuration updates are now defined 


M/ - Mi 

A t 


E - E 3i 

i , . — ” - 61 — 77 

Mm 1 ^ ” M3m* 


(3.7) 


Introducing eqs. (3.6) and (3.7) into the approximate expression for Hamilton’s Principle, eq. (3.3), then 
leads to 


1 j { 2 («». + C^™) • [(«/ - &) + COSs/ - £*)] + (v> - v-)} <Kdn 

= J J {% (*/ + <£»/) • (M/ + CC 3/ ) - f (Mi + c4i) ■ (Mi + CC,i) + (Vf - V)} dean 

= f f [( Kf - iQ + (Vf - Vi)] dCdn = 0 . (3.8) 

Jn Jh 


This result clearly implies the conservation of the total mechanical energy of the system within a step. 
In summary, the approximate form of Hamilton’s Principle given by eq. (3.3) leads to a discrete energy 
conservation statement, eq. (3.8), when the configuration updates are chosen according to eqs. (3.7), and the 
average stress according to eq. (3.5). 

Integrating through the thickness of the shell leads to 



r k f - hi 

At. 


(Vim, I + N. 2771 , 2 ) 

+ (£ 3 / — C.3i) 


At 


~ (Mj m> l + 


dn = 0 . 


(3.9) 


Tn this expression, the in-plane forces are N_ am — ( E m N* na + H m M ^ n ) / the out-of-plane forces V 3m = 
E m N ’ a , and the bending moments M_ nm = ( E m M * aa ) / . The discrete governing equations of motion 
for shells are then 


hf ~ hi 
At 


-(v, 


+ Kim, 2 ) ~ P m ' 


(3.10) 


Qr. 


kill 

At 


Qm (Mlm.l + M.2m,2 h^3m) ~ 2„,i 


(3.11) 


where p are the externally applied loads, and q * the externally applied moments measured in the local 
system. The finite change in director orientation E^f — E 3i was expressed in terms of the two parameter 
incremental rotation vector, see B. 
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Invariance of the system Hamiltonian under spatial translations arid rotations implies the conservation 
of the linear and angular momenta. Although discrete preservation of momenta is less crucial than discrete 
preservation of energy, it is interesting to note that eqs. (3.10) and (3.11) also imply the discrete preservation 
of this invariant. At first, eqs. (3.10) are projected onto the test functions k{Lq + lim) anr * ec l s - (3.1 1) onto 
the test, functions £E 3m , where 7r is an arbitrary vector. Next, integration over the shell reference surface 
yields 


J a {itfeo + M n 


h f ~h t 

At 


~ (N lm ,l + —'2m , 2 ) 
+ 




At 




dfi=0. (3.12) 


Straightforward algebraic manipulations then lead to 


^ • J [(£0 + M/) k f - (£0 + Mi) hi + k„Mj - Ui) 

+ E f g f - E tli + gJE f - !,)] dfi = 0, (3.13) 

where the following result was used 

—1 m — 1 m + ^ 2 m— 2 m h £ 3 h E3 t 2m^i-2m h ^3m^-3 m “ (3-14) 


Inserting the configuration updates, eqs. (3.7), into eq. (3.13) then yields 

■ j [(£ 0 + M/) kf - (to + Mi) ki + E y g f - E 3 + (k ni U,n + 5 m ^3 m )] Aft = 0. (3.15) 

It is easily verified that h m u rn + g m E^ m ~ 0. Hence, since 7r is arbitrary, eq. (3.15) implies the discrete 
conservation of the total angular momentum, J^[(r 0 + u) h + E_ gj dfi. Finally, projecting eqs. (3.10) onto 
the test functions 7r and eqs. (3.11) onto the null test functions gives the discrete conservation of the total 
linear momentum J^h dfi. 

It is important to note that any spatial discretization of the discrete equations of motion will inherit 
the discrete energy and momentum conservation statements just proved, when the configuration updates are 
chosen according to eqs. (3.7), and the average stress according to eq. (3.5). 


4. Energy Decaying Scheme. As discussed in the introduction, energy preservation, per is not 
sufficient to yield robust time integration schemes. High frequency numerical dissipation must be added as 
an inherent feature of the scheme. Such a scheme will now be constructed for the shell equations of motion 
using the EP scheme as a basic building block. 

First,, an additional state is introduced at time tj = lim £ _*o(£i + c), and the subscript (-)j is used to 
denote quantities at this time. The following averages are now' defined 

(•), = ![(■)/ + (■),•]; (-)/. = ![(); + (•)*]• ( 4 - 1 ) 

The ED scheme proceeds from the initial to the final time by means of two coupled steps: one step from t t 
to tf, the other from to tj. The time-discrete equations of dynamic equilibrium are 




T- g f ~ g , T 


(4.2) 
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(4.3) 


+ 3 [(^1.9,1 + ^29,2) - (i^ip.i + —29,2)] - p h \ 

Qi, + 3 (Ml 9 ,1 + M2 9 ,2 - — 3 9 ) “ Qh (Mip.l + M-2p,2 ~ — 3p ) ] = 

The configuration update relationships are given as 

Uj = u t + Af («, + n^/2, Uj = Mi - At [«/ - Mi - «(«j - Mi)] /6; 

£ 3 / = C.(, + A/ [E^f + Ej j ) / 2 , E^j = E 3i — At ~ ® 3 > — a {E-ij ~ 2£ii)] /®> 

where a is a tuning parameter that controls the amount of numerical dissipation provided by the scheme 
while the forces N_ np and moments M ap are given by 

K np =K nh +a (N aj - N ni )/ 2; M np = M nh + o (M aj - A£„,)/2. 


(4.4) 


(4.5) 


Using developments similar to those exposed for the EP scheme, it can be easily shown that the proposed 
discrete equations imply 


(Kf + V f ) - ( K, + Vi) +nc 2 =0. 


(4.6) 


c 2 is a positive quantity given by 



1 

2 

1 

2 


m || M II • It M II +2.9* || M II • II £3 II +T m || 4 II • II E 3 ||] dft 

IUII C* || e || dH > 0 , 


(4.7) 


where || • ||= (*)j - (•)* is the jump between ti and tj. This result implies the decay of the total mechanical 
energy over one step of the algorithm, (Kf + Vf) < (K t + U). The parameter a clearly controls the amount 
of energy that is dissipated within the step. Two such parameters could be used, controlling the amount of 
dissipated kinetic and strain energies, respectively, but this level of complexity does not seem to be necessary. 
The property of preservation of momentum observed in the EP case is lost in the ED algorithm. 

If the above ED scheme is applied to a single degree of freedom linear oscillator, the asymptotic value 
of the spectral radius of the amplification matrix, p 0 OJ is found to be p^ = (1 — a)/( 1 + a). For a = 1, 
p oc — 0, and asymptotic annihilation is achieved. If a = 0, p^ — 1, and in view of eq. (4.6), energy is 
exactly preserved. Hence, the ED scheme is in fact a family of schemes with a single tuning parameter, 
o, that controls the amount of high frequency numerical dissipation; both asymptotic annihilation or exact 
energy preservation can be achieved with the same scheme by using a — 1 or 0, respectively. 


5. Numerical Examples. All the examples described in this section will be treated with the proposed 
ED family of schemes corresponding to values of the tuning parameter a € [0, 1]. Although any value of a 
within this range can be used, the examples described here will contrast the two extreme choices. For a — 1 
(Poo — 0), asymptotic annihilation is obtained, and this will be called the ED scheme. On the other hand, 
for a — 0 (poo = 1), exact energy preservation is achieved, and this will be called the EP scheme. 

5.1. Clamped Half-Cylinder under Point Load. Consider the half-cylinder of radius B = 1.2 m 
and width b = 2 m depicted in fig. 5.1. The shell has a thickness t = 6 mm, is build-in along edge BC 
and free along the other. The structure is made of aluminum; Young’s modulus E = 73 GPa, Poisson’s 
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Fig. 5.1. Configuration of the clamped half- cylinder. 
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1 U 



t = 4.8 [sec] 


Fig. 5.2. Configuration of the system at various instants in time. 


ratio v — 0.30 and density p — 2700 kg/m 3 . At. point D, the shell is subjected to a concentrated load 
P_ = — Po(tf)(i, +i 2 + i 3 ). The magnitude of the load is 


Po(t) = 


P (1 — cos2nt/T)/2 t < T, 
0 t > 7 1 , 


( 5 . 1 ) 


where P = 0.1 kN and T = 2.0 s. The shell was modeled by a regular 8x4 mesh of quadratic elements. All 
simulations were run with a time step At — 5.0 1 0 03 s, for a total time of 6 s. 

Under the effect of the applied loads, the shell bends predominantly in the vertical direction, its direction 
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Fic. 5 . 3 . Displacement components at point D. U\: A; U-2 * °l f/3: □ . ED scheme: solid line; EP scheme: dashed line. 



Fig. 5.4. Tim,e history of bending moments at point M. M 1 1 : A; M\ 2- o; M22 • D. PD scheme: solid line ; EP scheme: 
dashed line. 

of least bending stiffness, as illustrated in fig. 5.2 that shows the configuration of the system at various instants 
in time. The three components of displacements at point D are shown in fig. 5.3; vertical displacements 
of up to 0.6 m are observed. The bending and twisting moments, Mn, M 22 , and M 12 , respectively, at 
point M are shown in fig. 5.4. Note the significant transverse and twisting moments associated with the 
three-dimensional motion of the shell. The components of in- plane and transverse shearing forces are shown 
in fig. 5.5 and 5.6, respectively. 

Next, the same problem was simulated with = 1, i.e. with no high frequency dissipation. The 
corresponding results are shown in figs. 5.3 to 5.6. Displacement and moment results are found to be in 
excellent agreement. At the scale of the figures, they are, in fact, indistinguishable. For the period of time 
2 < t < 6 s, the system is not. subjected to any loading, and the total mechanical energy of the system should 
remain constant. For the EP scheme, the energy is indeed preserved, as expected; for the ED scheme, 0.3% 
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Fig. 5.5. Time history of in-plane forces at point M. Fn ■' top /ignrp; Fi 2 : midci/e figure; F22* fco^tom figure. ED scheme: 
solid line; EP scheme: dashed line . 




Fig. 5.6. Time history of transverse shear forces at point M . F\ 3 ; top figure; F 23 . feoftom /igure. FT sc/temc: safari line; 
FP scheme: dashed line. 

of the energy is numerically dissipated in this period. It could be concluded that the EP and ED solutions 
are nearly identical, and that numerical dissipation is not necessary. However, the ED and EP scheme 
predictions for the in-plane and transverse shearing forces, shown in fig. 5.5 and 5.6. are markedly different.. 
EP predictions for force components F 12i and F\% show high frequency oscillations that are absent in 
the corresponding ED predictions. A simulation using the ED scheme with a time step At = 1.0 10“ 03 s 
showed that the ED predictions are converged. A simulation using the EP scheme and the same smaller time 
step yielded results with increased high frequency oscillations for the force predictions. It should be noted 
that the dynamic response of this simple system is very smooth; yet even here, high frequency numerical 
dissipation appears to be necessary to obtain a smooth, converged solution. 
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FlO. 5.8. Time history of the total mechanical energy (top figure) and relative energy loss (bottom figure). 
(EDS): solid line ; Plate model (EPS): dashed line; Beam model (EDS): dashed-dot line. 


Plate model 


5.2. Dynamic Response of a Plate with Edge Beams. Consider the rectangular plate of length 
L — 2 m, width b = 0.05 m and thickness t — 2 mm as depicted in fig. 5.7. Two circular beams of radius 
r — 2 mm are attached at the plate edges. A third circular beam is located at the center of the plate. All 
components are made of aluminum; Young’s modulus E — 73 GPa, density p — 2700 kg/m 3 . The total mass 
of the edge beams is 10 kg each, and that of the central beam is 1 kg. The plate is subjected to uniformly 
distributed loads F m and F/ along FE and CB, respectively. The components of these loads along the i } 
and z 3 axes are F 1/n — 40 N/m and F 3m = 80 N/m, respectively, and Fu = -20 N/m and F 3 / — -60 N/m, 
respectively. The common time history of each loading component is 


m) - 


Fi (1 — cos2nt/T)/2 t < F, 
0 t > T. 


(5.2) 
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FiO. 5.9. Trajectory of the plate mid-span point. Plate model (EDS): solid line; Plate model (EPS): dashed line; Ream 
model (EDS): dashed-dot, line. 



Fro. 5.10. Time history of the quarter-span axial force. Plate model (EDS): solid line; Plate model (EPS): dashed line; 
Ream model (EDS): dashed-dot line. 


where T = 3.0 s. The plate is discretized with 10 quadratic plate elements along its length. A constant time 
step Ai = 6 10~° 3 s was used for all simulations. 

At time t > 3 s the applied load vanishes, and the system total mechanical energy should remain 
constant. The evolution of the energy of the system for three different cases is shown in fig. 5.8: the plate 
model using the ED scheme, the same plate model using the EP scheme, and a simplified model of the system 
using beam elements and the ED scheme. For times t > 3 s, the total energy remains a constant for the EP 
scheme and nearly constant for the ED schemes. The relative energy loss is also presented in the figure: 0.6% 
of the energy was dissipated by the ED scheme in the period t £ [3,20] s. This figure clearly demonstrates 
the non-increasing property of the energy evolution for the proposed ED schemes, as opposed the constant 
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Fig. 5.11. Time history of the quarter-span transverse shear force. Plate model (EDS): solid line ; Plate model (EPS): 
dashed line; Beam model (EDS): dashed-dot. line. 

energy predicted in EP simulations. The trajectory of the plate mid-span point is shown in fig. 5.9: good 
correlation is observed between the predictions of the three models. The beam model is slightly off due to the 
inherent simplifying assumptions. The behavior of the quarter-span axial force and transverse shear force 
are shown in fig. 5.10 and 5.11, respectively. The poor predictions of the EP schemes are obvious in these 
two plots. The history of axial force presents violent oscillations with amplitudes an order of magnitude 
larger than those observed for the ED scheme. The history of the transverse shear force predicted by the 
EP scheme quickly diverges from the ED predictions for both beam and plate models. To ascertain the 
accuracy of the ED predictions, a convergence study was performed. Nearly identical results were found 
with smaller time step sizes At = 3.0 and 1.0 10 s, or when using the time adaptivity procedure. On 
the other hand, oscillations of increasing amplitude were found as the time step size is reduced in the EP 
scheme. Furthermore, the time adaptivity procedure failed to yield any results because the time step size 
was driven to unreasonably small values, At = lO -07 s, as the procedure tries to cope with increasingly 
violent oscillations. 
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Fig. 5.12. Configuration of the cruciform problem. 



Fig. 5.13. Configuration of the cruciform at times t — 0.062 and 0.093 s. 


5.3. Dynamic Response of a Cruciform. Consider a cruciform consisting of four thin panels ( Panels 
A, B, C ' and D) connected to a central beam, as depicted in fig. 5.12. Each panel is of thickness t — 4 mm, 
length L = 1.2 m, and width 6 = 0.1 m. The central beam has a square cross-section of width a = 8 mm. 
A mass M = 12 kg is attached at the tip of the central beam at point T. Panels and beam are simply 
supported at the root of the cruciform. A concentrated load P(t) is applied at point T. The load acts in 
the plane defined by axes i 2 and and makes a 30 degree angle with axis All components are made of 
aluminum with properties given in the previous example. The time history of the applied load is 


p(*) = 


Po (1 — cos27rt/T)/2 t < T, 
0 t > T, 


(5.3) 


where Pq — 1.2 kN and T = 0.1 s. 

As the applied load increases, in-plane stresses in the panels rapidly increase and buckling takes place 
in those panels subjected to compression, as can be observed in fig. 5.13 that depicts the configuration of 
the cruciform at two instants in time. The trajectory of point T projected onto plane i 2 , i 3 is shown in 
fig. 5.14. For reference, the corresponding trajectory of a beam with cross-sectional properties equivalent to 
those of the cruciform is also presented. Of course, the equivalent beam model is much stiffer since it does 
not allow buckling to take place. Furthermore, the motion remains confined to the plane defined by axis i } 
and the line of action of the applied load. When each panel is modeled individually, the stiffness of system 
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DISPLACEMENT U 2 [m] 


Fig. 5.14. Trajectory of point T of the cruciform projected in the plane i 2 , Solid line: shell model; dashed line: beam 
m.odel. 



FlG. 5.15. Total mechanical energy of the system. Solid line: shell model; dashed line : beam model. 

varies both spatially and temporally, giving rise to the more complex motion shown in fig. 5.14. The total 
mechanical energy of the system is shown in fig. 5.15. From time t — 0.1 to 0.2 s, the system is free and 
its total mechanical energy should remain constant. Due the dissipative nature of the integration scheme, a 
small amount of energy is dissipated over that, period of time: 2.7% of the energy was dissipated over the 
2435 time step period. 

The root, shear force arid quarter-span bending moment in Panels A and C are shown in fig. 5.16 
and 5.17, respectively. Each panel undergoes alternating phases of tensile and compressive loading. During 
the compressive phases, buckling takes place, and large shear forces and bending moments are observed in 
contrast with the tensile phases during which these quantities remain much smaller. 
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t = 0.06 [sec] * t = 0.12 [sec] 



L t = 0.18 [sec] t = 0.24 [sec] 


Fig. 5.19. System configurations at various time instants during the simulation. 


5.4. Snap-Through of a Cylindrical Shell. The snap-through behavior of a cylindrical shell under 
a concentrated load was investigated in ref [22]. The shell consists of a 60 degree sector of a cylinder of 
height h — 5 m, radius R = 5 m and thickness t — 0.1 m, as shown in fig. 5.18. Material properties arc: 
Young’s modulus E = 210 GPa, Poisson ratio v — 0.25 and density p = 10 4 kg/m 3 . The two straight edges 
of the shell are simply supported, while the two curved edges are free. 

A concentrated force F is applied at the shell’s apex. This force linearly increases from 0 to 5 10 7 N in 
0.2 s, then is held constant at that value. The simulation ends at time t — 0.3 s. Due to the symmetry of 
the problem, a quarter shell only is modeled; a regular 4x4 mesh of quadratic elements was used. The time 
step size was selected as At — 10“ 3 s. 

As the load increases, the shell apex displacement increases, then suddenly, snap-through takes place 
and curvature reverses. Curvature reversal initiates in the region of the applied load, then quickly propagates 
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TIME [sec] 

Fig 5.20. Time history of the plate center vertical displacement. ED schem.e: solid line; EP scheme: dashed line . 


x 10* 



TIME [sec] 

Fig. 5.21. Time history of the plate center forces for the ED and EP schemes. Force Fi : solid line; F 3 : dashed line. For 
clarity the EP results are shifted downwards 3 10 8 N. 

throughout the entire structure, which undergoes subsequent violent oscillations. Snapshots of the system 
at various instants in time are given in fig. 5.19. The vertical displacements of the point of application of the 
load computed with the ED scheme is shown in fig. 5.20. Note the gradual increase of the shell deflection, 
until collapse at buckling and the resulting vibratory response in the inverted configuration. 

Ref. [22] presents simulations of this problem using various schemes: the generalized-a [19] and the 
CEMA [22] schemes. The former scheme features high frequency numerical dissipation and linear stability 
properties, while the latter adds to the generalized-a method a constraint on the total mechanical energy of 
the system. CEMA is therefore both energy preserving and high frequency dissipative. The results presented 
in fig. 5.20 are in close agreement with those obtained with the generalized-a scheme, but quite different 


19 


1 i “ 1 



I I _ 1 _ I — 1 1 

0 005 0.1 0.15 02 0.25 

TIME [sec] 


Fig. 5.22. Time history of the plate center vertical velocity. ED scheme: solid line ; EP scheme: dashed line , shifted 
downwards 3 10 2 m/s for clarity. 


from those predicted by CEMA in the post buckling regime. It is important to realize that the higher 
modes are only an artifact of the discretization process, and should therefore be removed from the computed 
response. A standard scheme like the generalized-a method accomplishes this goal through the characteristic 
low-pass shape of its spectral radius; however, there is no guarantee that energy will not be allowed to grow 
within one step for nonlinear problems. In contrast, CEMA enforces the exact conservation of energy in the 
nonlinear regime, but at the same time inherits high frequency dissipation from the underlying generalized-a 
algorithm. Consequently, an artificial mechanism for transferring energy from the higher (artificial) modes 
to the lower modes is created that drives the response to an erroneous solution. In contrast, the proposed 
ED scheme achieves both nonlinear stability and high frequency dissipation. 

Next, the same problem w r as simulated with — 1, i.e. with no high frequency dissipation. In this 
case, two refinements in time step size were required to successfully complete the simulation, one at time 
t - 0.1665 s {At — 5 10~ 4 s), the other at time t = 0.2142 s {At = 2.5 10' 4 s). Deflections predicted 
by the EP and ED schemes, shown in fig. 5.20, are in good agreement during the initial snap-through 
phase, but become increasingly different during the subsequent oscillations. The force and velocity fields 
are markedly different. The plate center forces for both EP and ED schemes are shown in fig. 5.21. The 
forces predicted by the EP scheme present violent oscillations of amplitude up to an order of magnitude 
larger than those predicted by the ED scheme. These violent oscillations hamper the convergence of the 
Newton process at each time step, leading to the need for smaller time steps. The same observations can 
be made about fig. 5.22 which compares the plate center vertical velocity. Violent oscillations are initiated 
at snap-through and the strict preservation of energy implied by the EP scheme prevents any subsequent 
decay of these vibrations. Since vibratory stresses are a great importance to designers, it is essential to 
assess the ability of new integration schemes to reliably predict these quantities. It is unfortunate that many 
scientific publications about geometric integration only present responses for preserved quantities such as 
total mechanical energy or momentum. The above plots demonstrate that while EP scheme might perform 
very well for the prediction of total energy, momentum, or even displacement fields, they are unable to 
reliably predict other important fields such as velocities and internal stresses. Consequently, such schemes 
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are of little values in real life applications. 

6. Conclusions. Tn this work, a. robust algorithm for the dynamic analysis of geometrically exact shell 
structures was presented. The method is geometry-based, i.e. it incorporates knowledge about specific 
qualitative features of the underlying partial differential equations. However, departing from the classical 
approaches based on strict preservation of energy, the method presented here allows the system to drift away 
from the level set of constant energy in a controlled and tunable manner. 

This feature achieves two goals. First, a bound is placed on the total mechanical energy of the discrete 
system, leading to the concept of nonlinear unconditional stability; this stability criterion is stronger than 
that obtained through the classical analysis of numerical schemes. The resulting numerical procedure is 
endowed with superior robustness, an important feature when dealing with complex engineering problems. 
Second, the monotonic energy drift is associated with numerical dissipation of the high frequency modes. 
This tunable dissipation makes the algorithm stiffly accurate, and avoids the build up of energy in the higher 
modes that are an artifact of the spatial discretization process. 

The proposed scheme can deal with general shell structures and is not tied to a specific spatial dis- 
cretization of the governing partial differential equations. Kinematic nonlinearities are treated in a rigorous 
manner, and material nonlinearities can be handled when the constitutive law T s stem from the existence of a 
strain energy density function. The efficiency and robustness of the proposed approach w 7 ere demonstrated 
with specific numerical examples. 


REFERENCES 

[1] K. Bathe AND E. Dvorkin, A four-node plate bending element based on mindlin/reissner plate theory 

and a mixed interpolation , International Journal for Numerical Methods in Engineering, 21 (1985), 
pp. 367-383. 

[2] , A formulation of general shell elements - the use mixed interpolation of tensorial components , 

International Journal for Numerical Methods in Engineering, 22 (1986), pp. 697-722. 

[3] O. BauchaU, Computational schemes for flexible , nonlinear multi-body systems , Multibody System 

Dynamics, 2 (1998), pp. 169 225. 

[4] j On the modeling of friction and rolling in flexible multi-body systems , Mult ibody System Dynam- 

ics, 3 (1999), pp. 209-239. 

[5] , Analysis of flexible multi-body systems with intermittent contacts , Multibody System Dynamics, 

4 (2000), pp. 23 54. 

[6] O. BAUCHAU AND C. Bottasso, On the design of energy preserving and decaying schemes for flexible, 

nonlinear multi-body systems, Computer Methods in Applied Mechanics and Engineering, 169 (1999), 
pp. 61-79. 

[7] 5 Robust integration schem.es for flexible multibody systems , Computer Methods in Applied Me- 

chanics and Engineering, (2000). 

[8] O. Bauchau, J. Choi, and C. Bottasso, On the modeling of shells in multibody dynamics , Multibody 

System Dynamics, (2000). 

[9] O. BAUCHAU and T. Joo, Computational schemes for nonlinear elasto- dynamics, International Journal 

for Numerical Methods in Engineering, 45 (1999), pp. 693- 719. 

[10] O. Bauchau AND N. Theron, Energy decaying scheme for non-linear beam models , Computer Methods 
in Applied Mechanics and Engineering, 134 (1996), pp. 37-56. 


21 


[ 11 ] 


, Energy decaying scheme s for nonlinear elastic multi-body systems. Computers and Structures, 59 

(1996), pp. 317 331. 

[12] M. Borri, C. Bottasso, and L. Trainf.m.i, A novel momentum-preserving energy -decaying algo- 

rithm for finite element, mnltibody procedures , in Proceedings of Computational Aspects of Nonlinear 
Structural Systems with Large Rigid Body Motion, NATO Advanced Research Workshop, Pultusk, 
Poland, July 2-7, 2000. 

[13] M. Borri, L. Traineu.i, and C. Bottasso, On representations and parameterizations of motion, 

Multibody Systems Dynamics, 4 (2000), pp. 129 193. 

[14] C. Bottasso and M. Borri, Energy preserving/decaying schemes for non-linear beam dynamics using 

the helicoidal approximation , Computer Methods in Applied Mechanics and Engineering, 143 (1997), 
pp. 393 415. 

[15] , Integrating finite rotations , Computer Methods in Applied Mechanics and Engineering, 164 

(1998), pp. 307 331. 

[16] C. Bottasso, M. Borrt, and L. Trajnellt, Integration of elastic multibody systems using invariant 

preserving /dissipating schemes. Part I formulation & Part II num.eri.cal schemes arid applications. , 
Computer Methods in Applied Mechanics and Engineering, (2000). 

[17] M. BuOALEM AND K. BATHE, Higher-order mite general shell elements. International Journal for 

Numerical Methods in Engineering, 36 (1993), pp. 3729 3754. 

[18] C. Bunn and A. Isert.es, Geometric integration: Numerical solution of differential equations on m.an- 

j folds, Philosophical Transactions of the Royal Society of London Series A-Mathematical Physical 
and Engineering Sciences, 357 (1999), pp. 945-956. 

[19] J. Chung and G. Hulrf.rt, A time integration algorithm for structural dynamics with improved 

numerical dissipation: The generalized-a method , Journal of Applied Mechanics, 122 (1995), pp. 254 
266. 

[20] H. HlLBER, T. HUGHES, AND R. TAYLOR, Improved numerical dissipation for time integration algo- 

rithms in structural dynamics , Earthquake Engineering and Structural Dynamics, 5 (1977), pp. 282- 
292. 

[21] T. Kane and D. Levtnson, Dynamics : Theory and Applications , McGraw-Hill, Inc., New York, 1985. 

[22] D. KlJHL and R. E., Constraint energy momentum algorithm and its application to non-linear dynamics 

of shells, Computer Methods in Applied Mechanics and Engineering, 136 (1996), pp. 293 315. 

[23] J. SlMO AND N. Tarnow, The discrete energy-momentum method, conserving algorithms for nonlinear 

dynamics , ZAMP, 43 (1992), pp. 757-792. 

[24] , A . new energy and momentum conserving algorithm, for the nonlinear dynamics of shells , Inter- 

national Journal for Numerical Methods in Engineering, 37 (1994), pp. 2527-2549. 

[25] J. SlMO, N. Tarnow, and M. Doblare, Non-linear dynamics of three-dimensional rods: Exact energy 

and momentum conserving algorithms , International Journal of Numerical Methods in Engineering, 
38 (1995), pp. 1431 1473. 

[26] J. SlMO AND K. Wong, Unconditionally stable algorithms for rigid body dynamics that exactly preserve 

energy and momentum , International Journal for Numerical Methods in Engineering, 31 (1991), 
pp. 19-52. 


22 


Appendix A. Rodrigues Parameters. 

A common representation of finite rotations [21] is in terms of Rodrigues parameters r_ — 2k tan 0/2, 
where 0 is the magnitude of the finite rotation and k the components of the unit vector about, which it takes 
place. The following notation is introduced r 0 = cos 2 0/2 = 1 / (1 +r T r/4), and the finite rotation tensor 
R then writes 


R(r) = / + r 0 f+y rr. 


The following decomposition of the rotation tensor is extensively used in this work 


*=lr + - 2 


i + 


=HfH> Hf 


(A.l) 


(A. 2) 


Appendix B. Orientation of a Unit Director. 

Consider a unit vector called a director , that rotates to a final orientation For convenience, this 
director is considered to be the third unit vector of a triad S defined by ij , , £3, rotating to a triad S with 

orientation e x , e 2 , e 3 . The relationship between these two triads is e a = R i Q , where B is an orthogonal 
rotation tensor. If one solely focuses on the director, this rotation tensor is not uniquely defined, as any 
rotation about the director leaves its orientation unchanged. A virtual change in the director orientation is 


Je 3 = ej Sip , 

where Sip is the virtual rotation vector, Sip = 5RR T . 

The components of the virtual change in director orientation measured in S* become 


R r Se 3 - R r 7 % £0 = i£ R r Sip = IJSip* = 


-Sip* 

Sip; 

0 


(B-l) 


(B.2) 


where Sip* are the components of the virtual rotation vector in S * . This relationship clearly demonstrates 
that arbitrary values of Sip 3 , corresponding to virtual rotations of the director about its own orientation, 
will not affect virtual changes in the director orientation, and hence, setting Sip 3 = 0 is a valid choice. The 
following notation is adopted 


Sip* = i x Sn\ + uSaX — b Sa*\ b = [i.j , I2] - (B-3) 

So* is a 2 x 1, “two parameter” virtual rotation vector. It follows that £0 — R Sip* — Rb 5a * , and hence 

Se n = b 5a* . (B.4) 


If Rodrigues parameters are used to parameterize 7?, an equivalent expression can be obtained for finite 
changes in director orientation with the help of eq. (A. 2) 




i = Rm *1 i s‘ = Qm L* =b S', 


(B.5) 


where r* are the Rodrigues parameter measured in S * , and s* the corresponding “two parameter” incremental 
rotation vector. 
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